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Abstract 

We have performed a high statistics simulation of the 0(4) model on a three- 
dimensional lattice of linear extension L = 120 for small external fields H. Using 
the maximum entropy method we analyze the longitudinal and transverse plane spin 
correlation functions for T < T c and T > T c . In the transverse case we find for all 
T and H a single sharp peak in the spectral function, whose position defines the 
transverse mass rriT, the correlator is that of a free particle with mass ttit- In the 
longitudinal case we find in the very high temperature region also a single sharp peak 
in the spectrum. On approaching the critical point from above the peak broadens 
somewhat and at T c its position is at 2m<r for all our if— values. Below T c we 
find still a significant peak at u = 2rriT and at higher uo— values a continuum of 
states with several smaller peaks with decreasing heights. This finding is in accord 
with a relation of Patashinskii and Pokrovskii between the longitudinal and the 
transverse correlation functions. We test this relation and its range of applicability 
in the following. As a by-product we calculate critical exponents and amplitudes 
and confirm our former results. 
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1 Introduction 



0(N) spin models with N > 1 possess two types of two-point correlation functions, 
one for the transverse and and one for the longitudinal spin component, defined 

— * 

relative to the direction of the external field H. Correspondingly, there exist also 
two distinct susceptibilities and correlation lengths. The difference and interplay 
between the two sets of observables is of interest not only for the critical behaviour 
near to T c , the critical temperature, but as well for the study of the singularities 
induced by the existence of the massless Goldstone modes in these O(N) spin models 
for dimension 2 < d < 4 and all T < T c [1, 2]. 

In a previous paper [3] we studied the correlation lengths which govern the 
exponential decay of the two correlation functions in the three-dimensional 0(4)- 
model on lattices of linear extensions L = 48 — 120. We were able to confirm the 
predicted singular Goldstone behaviour [4, 5] of the transverse correlation length 
near the coexistence line, T < T c , H = and to determine the scaling function of 
£r. For the longitudinal correlation length the situation was different. In the high 
temperature region the scaling function could be calculated, in the low temperature 
phase, T < T c , we were however unable to reliably estimate £l, the data were not 
even scaling. This was ascribed to a spectrum of higher states which contribute to 
the longitudinal correlators below T c . In this paper we want to calculate the spectral 
functions of the two correlators and to actually find these states above the ground 
states rriT,L = V£t,l- Here, not only the spectrum of the longitudinal correlator is of 
interest, but as well the transverse spectrum. Spin- wave theory assumes, that long- 
wavelength transverse fluctuations dominate for small fields in thermal equilibrium 
below T c and that these fluctuations are describable by the Gaussian model or free 
field functional. This assumption can be tested readily by comparing the transverse 
correlator to the known Gaussian form. Of course, no higher state should then 
appear in the spectrum. 

The resulting dependence on H of the lowest states m T and m L at fixed T can 
be utilized to test critcal behaviour and the effects of massless Goldstone modes. At 
the critical point we had already found in Ref. [3] that m L = 2m T , like for 0(2) [6]. 
The equation between the two masses or correlation lengths follows from a relation 
between the transverse and longitudinal correlation functions [5, 4] which should be 
valid for small H in the whole low temperature phase. We shall discuss and test this 
relation explicitly for its range of applicability. In order to achieve these goals, we 
simulate the 0(4)-invariant nonlinear a-model with high statistics on a lattice with 
linear extension L = 120. We have chosen this specific model for several reasons: 
first because we want to clarify the open points raised in Ref. [3], second because 
in contrast to the corresponding 0(2)-model corrections to scaling are negligible 
here, and third because the model is of relevance for quantum chromodynamics 
(QCD) with two degenerate light-quark flavours at finite temperature, since its 
phase transition is supposed to belong to the same universality class as the chiral 
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transition of QCD [7, 8, 9]. 

The 0(iV)-model which we study here is defined by the Hamiltonian 

<x,y> x 

where x and y are nearest-neighbour sites on a three-dimensional hypercubic lattice 
with periodic boundary conditions, and (f>(x) is a iV-component spin vector at site 
x with length 1. It is convenient to give the following formulas still for general 
N, though we use in our simulations N — 4. The magnetization vector M is the 
expectation value of the lattice average of the local spins 

M = (0), with = I^0(f). (2) 

X 

Here, V = L 3 and L is the number of lattice points per direction, the lattice spacing 

— * 

a is fixed to 1. Due to the invariance of 7Y , the iJ-independent part of 7i, under 
0(iV)-rotations of <j>(x), the magnetization vector aligns with H 

M = Mn , with n = H/H . (3) 

— * 

It is therefore appropriate also to decompose the spin vectors <f>(x) into longitudinal 
(parallel to the magnetic field H) and transverse components 

0(f) = !l (f) ft + 4> L (x) . (4) 

The order parameter of the system, the magnetization M, is then the expectation 
value of the lattice average 0" of the longitudinal spin components 

M = ( 0" ) , and M x = ( ^ ) = . (5) 

Owing to the N components of M and H one is led to a whole matrix of suscepti- 
bilities 

Xab = d ^r b = d ~m =v (<^>-<^><^>) ' ^ 

with a,b — 1, . . . , N. The matrix can be expressed with only two quantities 

Xab = XLn a n b + xr^ab - n a n b ) , (7) 
the longitudinal and the transverse susceptibility. The total susceptibility is 

Xtot = Tr( X ) = Xl + Xt(N -1) = V( (0 2 > - (0) 2 ) . (8) 
From the last equation and 

XL = n X n = V ((0" 2 >-M 2 ) , (9) 



one finds, as expected, that the transverse susceptibility corresponds to the fluctu- 
ation per component of the lattice average (ft 1 - of the transverse spin components 

Xt = j^($ ±2 ). (10) 
The two susceptibilities can as well be obtained from the equations 

dM a M m\ 

The remainder of the paper is organized as follows. We start with a detailed 
discussion of the correlation functions and their properties, the formula connecting 
the two types of correlators and their spectral representations. In Section 3 we 
address the critical behaviour of the observables. In the following Section 4 we 
describe some details of our simulations and the used maximum entropy method 
(MEM). Then we analyze our data in the different T-regions: below T c , at T c and 
in the high temperature region; first separately for the transverse and longitudinal 
correlation functions, subsequently we investigate their interplay. We close with a 
summary and the conclusions. 



2 The correlation functions 
2.1 2-Point correlators 

As for the susceptibilities, there exists a matrix of 2-point connected correlation 
functions 

Gab(x U X 2 ) = (</>a(Xi)<f> b (x 2 ))c = (0a(^l )<Pb{.X 2 ) ) ~ (0 a (f l) ) (0 6 (f 2 ) ) • (12) 

Furthermore, because of translation invariance we have 

(0(f)) = (0) , and Gab(xi,x 2 ) = G ab {x 2 - x i) . (13) 
The correlation functions are related to the susceptibilities by 

x 

The matrix G must have the same form as the matrix x 

Gab = G L n a n b + G T (5 ab - n a n b ) , (15) 

leading to 

G L (x) = (0"(O)0ll(f))-(0ll) 2 , (16) 
Gt{x) = j^i^m^x)). (17) 



As mentioned already, the large distance behaviour of the correlation functions Gl,t 
is determined by the respective exponential correlation lenghts ^l,t 

G L , r (£)~exp(-|£|/£ iiT ) . (18) 

This is true for all temperatures T (the coupling J acts here as inverse temperature, 
that is J = 1/T) and external fields H, except on the coexistence line H = 0, T < T c 
and at the critical point, where the correlation lengths diverge. 



2.2 2-Plane correlators 

Instead of the 2-point correlation functions we shall actually measure 2-plane cor- 
relation functions. They are defined as the connected correlation functions of spin 
averages over planes. For example, the spin average over the (x, |/)-plane at position 
z is 

The general 2-plane connected correlation functions are then 

D ab (z) = J2G ab (x) = L 2 {S a (0)S b (z)) c . (20) 

Also the expectation values of the longitudinal and transverse components of S(z) 
are independent of z and equal to those of the corresponding lattice averages 

(Sll(z)> = (0">=M, and (S ± (z)) = (^)=0. (21) 

The respective longitudinal and transverse plane-correlation functions D LjT (z) are 

D L {z) = L 2 ((S\0)S\z))-M 2 ) , (22) 

D T (z) = JL-(^(Q)§^ z )) . (23) 

Here, z is the distance between the two planes. Instead of choosing the z-direction 
as normal to the plane one can as well take the x- or y-directions. Accordingly, we 
enhance the accuracy of the correlation function data by averaging over all three 
directions and all possible translations. The correlators are symmetric and periodic 
functions of the distance r between the planes, the factor L 2 on the right-hand sides 
of (22) and (23) ensures the relation 

L-l 

XLtT = J2D L Ar) ■ (24) 

r=0 

Like the point-correlation functions in Eq. (18) the plane correlators Dl,t(t) decay 
exponentially. 
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2.3 The Gaussian model 



We noted already in the introduction, that the Gaussian model is assumed to de- 
scribe the long- wavelength transverse fluctuations. Therefore, we shortly discuss, 
following Ref. [10], the aspects of the model which are relevant for our work, first 
in d dimensions, finite volume and the continuum limit. Here, 0(f) is a single 
component, real scalar field with the Hamiltonian 

(3Hg = J d d x[ ic s (V0) 2 + l -a s( t> 2 ] , (25) 

where c s is the so-called stiffness and a s is proportional to the square of the mass of 
the field. Let us introduce the Fourier transform 4>(k) of the field 

4>{k) = J d d xe- l%£ (t){x) , (26) 

m = ^x> ifc >(£). (27) 

k 

For periodic boundary conditions, the components of k are restricted to fcj = 27mj/L 

~ — * ~ — * ~ — * 

with rii = 0, ±1, . . . , and because is real we have (f)(k)* = <f)(—k). In terms of <p(k) 
the Hamiltonian is 

^ = ^^|0(£)| 2 (c s P + a s ) , (28) 



2V 

k 



~ — * 

that is the <f>(k) decouple, 



<0(£)> = O, and (|0(fc)| 2 ) = ^ ■ (29) 

CgPj ~\~ Ob g 

That leads to the results 

(0(f)) = 0, and (0(O)0(f)) = l^ f* , (30) 

v r c s k 2 + a s 

k 

The latter expectation value is thus directly the 2-point connected correlator of 
0(f). It is easy to evaluate the Gaussian correlator for d = 1, the case of the 2-plane 
correlation function 

D ( T ) = jX^ c s kl + a s = Lc~ s ^ kl+m 2 ' (31) 



Here, k n = 2nn/L and the mass is m = (a s /c s ) l l 2 . Using the well-known relation 

1 ^ e lknT e~ XT 

x + ik n ~ 1 - e' Lx ' 



1 piknT p XT 

1 _! = _1 (32) 



n 
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for x = ±m, we obtain 

i p -rriT I „-m(L-T) 

D(r) = --- • , f . (33) 

v ; 2mc s 1 - e- mi v ; 

In addition, we can identify a s with the inverse of the susceptibility 

L 

X = J D(r)dr = -i- = - . (34) 
o 



m 2 c s a s 



In the infinite volume limit, V — > oo, the Gaussian correlator is 



oo 



1 f d^k gikx _^ y y 

G ^ ~~ s J (2^ p + m 2 ' aIld G(A;) ~ ^ p + m 2 • 



-OO 



The integral is known for general d [10, 11] 

d/2-1 



G(f) = (i) _ ^-H2|). (36) 

where K v (z) is a modified Bessel function. We are interested here in the two cases 

e -m\x\ 
-mr 

D(r) = — , i=l. (38) 

The last equation can also be obtained from Eq. (33) for L — > oo. 



2.4 The relation between the transverse and longitudinal 
correlation functions 

The connection between the two correlation functions was first established in a paper 
by Patashinskii and Pokrovskii [5] . The main arguments leading to this relation are 
the following. The transverse susceptibility xt, Eq- obviously becomes infinite 
for all T < T c and H — > 0, due to the non-zero spontaneous magnetization M in 
the low temperature phase. The instability comes about because of infinitesimal 

— * 

long-wavelength transverse fluctuations which change the direction of but not its 
length 0. For a fluctuation 

5<p = - (0) , 50" = 0" - M , = $ L , (39) 

the length change is 

50 2 = 5<f) 2 = 2 - M 2 = (50") 2 + (50 x ) 2 + 2M50 11 . (40) 
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To lowest order in 50" length conserving fluctuations 

50 2 « (50 x ) 2 + 2M50 11 = , (41) 
entail longitudinal fluctuations of size 

r 2M v ; 

Patashinskii and Prokrovskii extend this "principle of conservation of the modulus" 
also to local fluctuations, that is to non-uniform spin-waves with long wavelengths 

The longitudinal correlation function Gl(x) can then be expressed as 

G L (x) = (0 11 (0)0" (x)) c = (50 11 (0)50" (£)) (44) 
1 '{ 50^(0) )\5$\x)) 2 ) . (45) 



4M 2 



c 



In the last expectation value we have explicitly noted that the connected value has 
to be taken. As already mentioned, the transverse fluctuations 50 ± (x) = 0- L (af) are 
assumed to have the Gaussian form of Eq. (25) 

(3H T = Jd 3 x[ \c s {V^f + ^-0 ±2 ] . (46) 

The stiffness c s depends on the temperature and is to lowest order independent of 

H . It measures the resistance of the system against non-uniform rotations of <f>(x). 

In the disordered phase, where the correlation length is finite, distant points can 

have different orientations and c s should be essentially 1. From Eq. (46) it is clear 

— * 

that the components of do not interact, since 

3 N / All \ 2 

^)'=EE| (47) 

i=l a =l v 1 7 

and the components 0f decouple. We may therefore use Wick's theorem to reduce 
the expectation value in Eq. (45) 

(0f(O)0f {x)) = 2<0^(O)0^)> 2 + <0f (O))(0f (£)> , (48) 

so that 

JV N 

= S <#(°)#(*)> 2 = 2M2 E ( WOW* ~ ^n b ) f , (49) 

a,b=l a, 6=1 
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where we have used Eq. (15). That finally leads to 

Gl(x) = ^ GUx) . (50) 

Let us briefly discuss the consequences of Eq. (50), the PP-relation, in the infinite 
volume limit. The transverse correlation function is then given by Eq. (37) with the 
mass m T and 

m T = = T7 — , (51) 

XTC S Mc s ' 

whereas the longitudinal correlation function is 

AT — 1 p -2mr\x\ 

G ^> = W (i^W ■ (52) 

The Fourier transform of Gl is 

GAk) = ^r- r } - \, arctan - - , (53) 
y ' 2M 2 4irc 2 s k 2m T ' v ; 

where k — \k\. It is instructive to compute the corresponding longitudinal suscepti- 
bility 

» = e *<*=°> = widb' (54) 

or 

to = l4i^F'^ 1/2 ' (55) 

The last equation implies the divergence of xl ~ H^ 1 / 2 for H — > and T < T c . 
That is exactly the singular Goldstone behaviour predicted in Refs. [12] and [1, 2] 
from field theory. Furthermore, we find from Eq. (51) that near the coexistence line 

m T ~ H 1/2 , or & ~ H' 1/2 , (56) 

and from Eq. (50), if Gt{x) is decaying exponentially 

6r = 2&. (57) 



In order to test all these relations with our data we need the lattice versions of the 
one-pole correlation functions, Eqs. (37) and (38), and the corresponding predictions 
from the PP-relation (50). On a lattice with periodic boundary conditions the 
Fourier transforms are 

= J? Y,^ s G{k) , G(k) = Xy^G(f) , (58) 

k $ 
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where x = (711,712,713) with n« = 0, . . . , L — 1 and k = 27r(mi, 7712, m%)jL with 
TOj = 0, ±1, . . . , ±(L/2 — l),L/2 for even L. Evidently, the 2-plane correlation 
function D T is proportional to the sum of exponentials 

D T ( T ) ~( e- mTT + e- mT(L " r) ) . (59) 

The normalization factor is determined by Eq. (24), so that 

/m,T\ e ~ mTT + e - m T(L-r) 
D T (r) = XT tanh (Tf) • . (60) 

For small itit the prefactor 

k = Xt tanh (^^j (61) 
coincides with that of Eq. (33). The Fourier transform of Dt is 

^ /,x 2sinh(mT) 

D T (k) = k ^ — / ; o~~FT\ • 62 

V ; 4sinh 2 (^) +4 sin 2 (|) V ; 

As in Ref. [13] we extend the form (62) to the three-dimensional Fourier transform 

n (T\ 2sinh(m T ) 

G T (k) = k , (63) 

4sinh 2 (^) + 4£sin 2 (|) 

i=i 

from which Gt{x) is obtained, using Eq. (58). Once we have determined m*r from 
MEM or a direct fit of the transverse 2-plane correlator data to Eq. (60), we can 
compare the longitudinal data to the following predictions from the PP-relation. 
Let us first compute the longitudinal susceptibility 

^ = (64) 

X 

From 

= ZsEE^&.o G T (k)G T (k') = L^G T (k)G T (-k) , (65) 

x k k' k 

where the prime on the sum indicates that there is no contribution from hi = n , we 
find 

» = WI5E'^(*) (66) 

k 



N - 1 
2/f 2 L 3 



s ^(f)E'U 2 (f) + E- 2 (!) 

c i=i \ / 



(67) 



9 



In the last equation we have used Eq. (61) and xt = M/H. Similarly, we can find 
the longitudinal 2-plane correlator from 



(68) 



x,y 



Here, we have 



Dl(z) 



or 



N-l 

2M 2 L 4 



e l ^ +k '^Y,'GT(k x ,k y ,k z )G T (k x ,k y ,k' z ) , (69) 



k z ,k' z 



k x ,ky 



fl 'W = ^^(^)E^E' ^(f ) + 5> 2 (|) 

k 3 ,k' 3 fci,fc2 L i=l 



X 



sinh 2 +^sin 2 ^ 

i=i ^ 



sin 



2 / ^3 



(70) 



2.5 The spectral functions 



We define the spectral function A(uS) of the 2-plane correlation function by the 
integral equation 



D(t) = / duuA(uu)K(uo,T) , 



(71) 



where the integral kernel is 



K(u,t) 



uit _|_ g— w(L— r) 

1 - e-^ L 



cosh(w(r - L/2)) 
sinh(cuL/2) 



(72) 



This choice of the kernel is motivated on one hand by the lattice form of the 2- 
plane correlators as obtained in the transfer matrix (TM) formalism [13] and second 
from the appearance of the same factor K{uj,t) in the continuum, Eq. (33), and 
lattice forms, Eq. (60), of the one-pole correlation function. Moreover, this kernel 
corresponds to the one known from finite temperature QCD [14, 15], where L is 
replaced by (3 = 1/T. Yet, we are free to redefine the kernel and the spectral function 
by a convenient factor q(ou) and its inverse such that D(t) remains invariant. In a 
paper by Aarts et al. [16] this freedom has been used to remove numerical instabilities 
in Bryan's MEM- algorithm [17] at small oo which are due to the singularity of the 
kernel at u = by changing the kernel with the factor (in our notation) 



q(u) = tuL/2 , K = qK , A — A/q . 



(73) 
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Here, as for the original kernel K(u>,t), the limits L — > oo and uj — > of .K" do 
not commute, because of the additional factor L in g. Unfortunately, not only the 
small but also the large ^-dependence of the kernel is changed by the choice (73). 
We arrive at a more adequate redefinition by first noting that for finite L in the 
continuum 



/ 



drK(u, r) = - , (74) 

UJ 



and on the lattice 



The choice 



L-l 

^K(^,r) = coth| . (75) 



r=0 



. uj/2 in the continuum 

q(uj) = <{ (76) 
tanh(a;/2) on the lattice , 

normalizes then K(u>, r) in r to 1, independent of L, and leads to the sum rule 

oo 

X = J dujA{uj) . (77) 



This approach has several advantages: it cures the numerical problems appearing in 
Bryan's algorithm, the limits of K for L — > oo and uj — > are finite and commute, 
and the large u;-behaviour remains unchanged on the lattice. From the measured x 
the sum rule gives us already some information on the size of A{uj) before we start 
to calculate it. We note here, that the bosonic kernel had been redefined earlier by 
Jarrel and Gubernatis [18] with a factor uj to ensure positive definite spectra. 

In applying the maximum entropy method a default model for the spectral func- 
tion is required which represents prior knowledge about A{uj) and accordingly A{uj) . 
Quite obviously the transverse default model is always 

A?(uj) = X TS(uJ-m T ) . (78) 

The longitudinal case is more involved. In the high temperature region a single pole 
ansatz as for the transverse spectral function seems to be adequate, since there for 
small H all spin components must show the same behaviour. For the low tempera- 
ture phase we may derive a model spectral function directly from the considerations 
in the previous subsection. The corresponding lattice version of the spectral function 
from Eqs. (69) and (70) is rather bulky. However, in the continuum limit we can 
get a useful formula for Dl(t) from the PP-relation in the same way as we obtained 
Eq. (33). The result is 



n ( \ N ~ l 1 \^ 



2M 2 c 2 s L 2 

&x fey 



z -m k T _|_ e -m k (L-r) 

2m* (1 - e~ m " L ) 



(79) 
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where m 2 , = k\ + A; 2 + m\ . This can be rewritten as 



D ( )- iV ~ 1 1 V 1 
MTj 2M 2 c 2 L 2 ^ (2m fc ) 2 

&x fey 



coth r ^ ^ (2mfc ' r) 



+ 



V 2 J m k 2sinh 2 (^) 



(80) 

which shows that Al consists of 5-function contributions at u = Im^ > 2m^ and at 
id — 0. The latter disappears exponentially with increasing L . The infinite volume 
expression is easily derived from Eq. (79) 



CO 



TV — 1 /• rf 2 A; e" 2mfcT iV — 1 /" , e- 2mfcT " 
DtAt) = / ^ — s- = — , o / am* 



2M 2 c 2 7 (2tt) 2 4m 2 16ttM 2 c 2 y m fc 
N-l 



-co mj 



16vrM 2 c 2 y w 



rfcu . (81) 



Since the kernel reduces for L — > oo to 



K(c,r) = |e— , (82) 



we obtain a continuous spectrum starting at u = 2rriT and 



AL{u) ~MPd- ^ • (83) 



The sumrule (77) leads again to the result (54) for xl ■ 



3 Critical behaviour and scaling functions 

In the thermodynamic limit (V — > oo) the observables show power law behaviour 
close to T c . It is described by critical amplitudes and exponents of the reduced 
temperature t — (T — T c )/T c or magnetic field H . The scaling laws at H = are 

for the magnetization 

M = B(-tf for t < , (84) 
the longitudinal susceptibility 

Xl = C + t^ for t > , (85) 

and since for H = 0, t > the correlation lengths coincide £ T = £l = £ (like the 
susceptibilities) 

f = f +r" for t > . (86) 
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On the critical line T = T c or t = we have for H > the scaling laws 

M = B C H 1/S or H = H M S , H = (B c y 5 , (87) 

XL = C c H l/5 - 1 with C c = B c /5 , (88) 
and for the correlation lengths 

Zl,t = Cl,tH~ Uc , v c = is/PS . (89) 

We assume the following hyperscaling relations among the critical exponents to be 
valid 

7 = P(6-1), dv = (3(1 + 5), 2-77 = 7/1/. (90) 

As a consequence only two critical exponents have to be known. The dependence 
of the observables on the temperature and the external field in the critical region is 
described by scaling functions 

M = hV'Mz) , XL = h^-'f^/Ho , (91) 
= h-^gf(z) , a = h-'-gfrz) . (92) 

Here, the scaling variable z is given by 

z = ih~ 1/A , A = PS , (93) 

where t = tB 1 ^ 13 , h = H/Hq are the normalized reduced temperature and field. The 
normalization is chosen such that 

/ G (0) = 1 , and fQ(z) x =_J-z)f> . (94) 

The scaling functions /g(-z) and f x (z) are universal and have been discussed in detail 
in Ref. [3]. The functions g^' (z) are universal except for a normalization factor. 
On the critical line t = or z = we find from (89) 

0€ ,T (O) = a, T (Br /fi , (95) 
and from (86) the asymptotic behaviour at z — > 00, that is for H — > 0, t > 

g^(z) =^ ?Wl*z- v . (96) 

The ratios of the amplitude for z — > 00 in (96) and the (0) are universal and 
therefore also the normalized scaling functions 

gf{z) = g L ^ T {z)/gf{Q) . (97) 

The stiffness c s in the critical region may be derived from Eq. (51), using ttit = 1/£t 
and the scaling functions 
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The leading term of c s in the high temperature region H — > 0, t > can be calculated 
in terms of critical amplitudes from Eq. (96) and the corresponding equation for 

/g(*)[3] 

fo( z ) 2 =oo C + H Q B & - l z^ = R x z^ . (99) 
We find the simple result 

Cs = (e+) 2 (c + ) _1 r"" • [i + o(r 2A H 2 )] . (ioo) 

Here, the higher terms are proportional to if 2 because £t is an even and M an odd 
function of H [3, 19]. Therefore, c s is given for small H essentially by the leading 
term and its value is close to one. In the low temperature phase both £ T and xt 
diverge because of the Goldstone effect. A corresponding analysis leads to the form 

c s ~ (-*)-"" ■ [1 + O ( (-t)- A/2 # 1/2 )] , (101) 

and here for small H the higher terms cannot be neglected, they even change the 
temperature dependence of c s . We note, that the stiffness p s defined by Privman et 
al. [20] is related to our definition by p s = M 2 c s . 



4 Numerical details 



All our simulations were done on three-dimensional lattices with periodic boundary 
conditions and linear extension L = 120. As in Ref. [3] we have used the Wolff single 
cluster algorithm [21], which was modified to include a non-zero magnetic field [22]. 
The coupling constant region which we have explored was 0.90 < J < 1.2, the 
magnetic field was varied from H = 0.0001 to H = 0.007. We performed 150-400 
cluster updates between two measurements, such that the integrated autocorrelation 
time of the magnetization was typically of the order of 1. Only at the smallest 
H— values Ti nt reached eventually 5. In general we made 100000 measurements at 
each fixed H and J. Besides M we measured the susceptibilities Xl,t and the 
correlation functions (^"(O)^'^)) and (S ± (0)S ± (r)) of the spin plane averages for 
the plane distances r = 0, 1, . . . , L/2. Before the start of the measurements 15000- 
20000 cluster updates were carried out to thermalize the system. In addition, up 
to 5000 of the first measurements were discarded as long as the average cluster 
size was still rising. In order to save storage we grouped the correlation function 
results in blocks of 500 measurements. The block averages led thus to 190-200 sets 
of correlation data for the transverse and the longitudinal spins respectively. In 
the latter case we completed the data to connected correlation functions using the 
magnetization averages of the corresponding block of updates. These data were 
finally used as input to our MEM-program. 

In the following we are using J c = T" 1 = 0.93590 from Ref. [23] as a value for the 
critical coupling. 
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4.1 The application of MEM 



As already indicated in Section 2.5 we determine the spectral functions with Bryan's 
MEM-algorithm [17] and the modified kernel K(uj, r) where the lattice version of 
q{oj) = tanh(u;/2) is used. The Shannon entropy term which contains the prior 
information is then 

oo 

S = J duo [A(uj) - A d (uj) - A(uo)\n(A(uj)/A D (uo))] , (102) 
o 

where A d (uj) represents the default model. The MEM-algorithm minimizes the 'free 
energy' 

F = X 2 /2 - aS , (103) 

or equivalently maximizes the conditional probability P(A\DA D ) ~ exp(— F). Here, 
X 2 /2 is the logarithm of the likelihood function for the data D (assuming a Gaussian 
distribution of the data) and a is a parameter balancing the relative importance of 
the data and prior knowledge. Bryan's algorithm calculates the spectral function at 
fixed a in a whole range of a-values and estimates each time a conditional probability 
P(a\DA D ), which is typically sharply peaked at some point a. The final result for 
A(uj) is then obtained by integrating over a with the conditional probability as 
weight. Further details of the MEM data analysis are described in Refs. [14] and 
[17]. Our MEM program is based on an earlier version which was written by Ines 
Wetzorke and was already used in Ref. [15]. 

Due to the large number of J and H values, where we simulated the 0(4)-model, 
namely 95, we obtained rather different correlators. In order to give an impression 
of our data we show the correlators Dl(t) for the longitudinal sector at J = 0.92, J c 
and J = 0.95 and all if -values in Figs. 1, 2 and 3, that is one plot each for T > T c , at 
T c and for T < T c . The correlation functions were calculated from all data belonging 
to the same T and H pair using the jackknife method. The results coincided within 
the errorbars with the averages from our MEM-input data sets. We observe at all 
couplings for increasing H an increasing descent of the correlators and accordingly 
a shorter range of r where the data are precise enough to be relevant for the MEM- 
analysis. Whereas in the high temperature region the logarithm of the correlator 
is essentially a straight line, that is we have a sharp single peak in the spectral 
function, the behaviour is gradually changing with T to a more complex spectral 
function in the low temperature domain. 

We proceed in the following way. Before we start with the MEM analysis of a 
transverse or longitudinal correlator we check whether it can be fitted with a one- 
pole ansatz, Eq. (60) (in the longitudinal case evidently with xl and m^). If the fit 
is not completely satisfying we start MEM by choosing a default model, eventually 
peaked at the previously found pole position. The suitable form of the default model 
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Figure 1: The correlator D l (t) at J = 0.92 for H = 0.0001 (circles), 0.0002 (tri- 
angles), 0.0005 (davidstars), 0.001 (triangles down), 0.002 (pluses), 0.003 (boxes), 
0.004 (filled triangles), 0.005 (diamonds) and 0.007 (filled boxes). 
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Figure 2: The correlator Dl(t) at J = J c for the same if- values as in Fig. 1 and in 
addition for H = 0.0003 (stars), 0.0007 (crosses) and 0.0015 (filled circles). 
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Figure 3: The correlator Dl(t) at J = 0.95 for the same if-values and with the 
same notation as in Fig. 1. 



as such is not the only important requirement for a sensible result of the analysis. 
The average size of the default model is also of significance. We achieve this by 
normalizing the default model such that 

X = J dcoA D (co) , (104) 
o 

where u m is the maximal value of u we use in MEM. If the integral in the last 
equation is definitely different from x then the form of the default model is irrelevant 
for the MEM result, we shall find only solutions for small a- values with a peak at a 
in the range of 0.5 to 10 and a spectral function consisting of isolated sharp peaks. 
Indeed, MEM performs then essentially only a x 2- fit) because the entropy term is 
suppressed by the small a. The choice of the cu-range is also of importance. Here we 
should keep in mind, that the spectral function at small uj determines the behaviour 
of the correlator at large r and vice versa the large u-region influences the correlator 
at small r. If the selected maximal value u m is too small MEM delivers nevertheless 
a spectral function which leads to a good correlator fit, however at the expense of 
a spurious sharp peak in the spectral function at uj m and a distortion of the peaks 
at lower u. Alternatively, one may omit the lowest r- values from the analysis, but 
then no information for large u is obtained. Therefore we just extend the tu-range 
till the peak at the end point disappears. We choose the a-range in such a way that 
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the conditional probability P(a\DA D ) at the lower and upper interval limit is about 
1% of the value at the maximum position a. Finally we check the MEM result for 
the spectral function by reconstructing the correlator and comparing the result with 
the correlator data. 



5 Results for the spectral functions 
5.1 The transverse spectrum 

As a first check we have tried a fit of our transverse correlator data with the Gaussian 
form, Eq. (60). The value for xt required in the equation, was either calculated from 
Eq. (11) or Eq. (24). Also, the r-range for the fit was varied to obtain sensible results. 
In all cases, that is for all our J and if -values, the transverse correlator data could 
be represented by the Gaussian one-pole form, Eq. (60). As examples we show in 
Figs. 4 and 5 our transverse correlator data for H = 0.0001. In Fig. 4 all data for 
T < T c are shown, in Fig. 5 the data at T c and for T > T c . Below T c the m^-values 
for H = 0.0001 are rather small so that the correlators are comparatively flat, the 
overall size, being proportional to xt — M/H, diminishes when T c is approached 




10 20 30 40 50 T 60 

Figure 4: The correlators D t (t) for H = 0.0001 below T c at J = 0.94 (circles), 0.95 
(triangles), 0.97 (boxes), 1.0 (diamonds) and 1.2 (stars). The lines are the respective 
fits with the Gaussian form, Eq. (60). 
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Figure 5: The correlators Dt(t) for H = 0.0001 at T c (filled circles) and above T c 
at J = 0.93 (triangles), 0.92 (boxes), 0.91 (diamonds) and 0.90 (stars). The lines 
are the respective fits with the Gaussian form, Eq. (60). 

from below. Above T c the average size diminishes further with increasing T, however 
the mr- value increases (the correlation length £t — 1/Wr decreases) such that the 
correlators become steeper. 

In Figs. 6 and 7 we show the results from the Gaussian fits for the transverse mass 
rriT below and above T c . The data points for the same J-values have been connected 
by straight lines to guide the eye. Since we expect in the low temperature region 
near the coexistence line that my ~ H 1 ^ 2 (see Eq. (56)), we have plotted vtit versus 
if 1//2 in Fig. 6. Indeed, we observe with decreasing temperature an increasing range 
in H 1 / 2 with a straight line behaviour. On the other hand, close to T c , here for 
J = 0.94, the curvature of rriTiH 1 ^ 2 ) grows slightly. Above T c we expect for H — ► 
a finite transverse correlation length and therefore that reaches a finite value. 
This expectation is borne out in Fig. 7, where we have plotted mx as a function of 
H. As in the low temperature range the curvature of itlt{H) is growing when T c is 
approached. The if-dependence of the transverse mass on the critical line T = T c 
results from Eq. (89) in 

m T = ^r 1 . H Vc . (105) 

In Fig. 8 we show the m^-data at T c as a function of H u °, where we used the value 
u c = 0.4024 obtained in Ref. [3]. The line through the data points is given by Eq. 
(105) with = 0.838 again from [3]. The data fully agree with the expectations. 
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Figure 6: The transverse mass my as a function of H 1 ^ 2 in the low temperature 
phase T < T c . The J- values are 0.94 (circles), 0.95 (diamonds), 0.97 (squares), 1.0 
(triangles) and 1.2 (davidstars). 
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Figure 7: The transverse mass m>T{H) for T > T c . The J-values are 0.93 (circles), 
0.92 (diamonds), 0.91 (squares) and 0.90 (triangles). 
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Figure 8: The transverse mass Tot (circles) versus H Uc on the critical line. The 
straight line is given by Eq. (105) with the parameters from Ref. [3]. 

A straight line fit to the data for In ttlt as a function of In H leads to the following 
results for the parameters 



confirming the previous numbers. 

5.2 The longitudinal spectrum 

In the high temperature phase the longitudinal and transverse correlators must coin- 
cide for H — > and therefore as well their spectral functions, that is the correlators 
will both have the same Gaussian form. With increasing H and/or decreasing tem- 
perature the two correlators will increasingly differ until the longitudinal mass ttil 
reaches 2m^, the threshold value of the longitudinal continuum in the low temper- 
ature phase. Indeed, that is what we observe in our analysis. 

At the two highest temperatures, corresponding to J = 0.90 and 0.91, we find that 
the longitudinal correlators for all our if- values are of the Gaussian form, Eq. (60). 
In Figs. 9 and 10 we show the results for the longitudinal mass in comparison 
to the respective transverse mass rrij- and 2rriT as a function of H. We see that 
at the lowest if- value the two masses coincide and that rriL becomes larger than 
tut with increasing external field. The relative difference increases with lowering 
temperature. In approaching T c further, at J = 0.92 and 0.93, we still find a single 
peak in the spectral function, however instead of a 5-function contribution, the 



£t = 0.839 ± 0.004 , and v c = 0.4020 ± 0.0007 , 



(106) 
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Figure 9: The longitudinal mass mi{H) for J = 0.90 (filled circles). For comparison 
ttlt and 2m*r are also shown. 
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Figure 10: The longitudinal mass ttil(H) for J = 0.91 (filled circles). For comparison 
rriT and 2irtT are also shown. 
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Figure 11: The longitudinal mass ttil(H) for J = 0.92 (filled circles). For comparison 
rriT and 2m*r are also shown. 
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Figure 12: The longitudinal mass m^H) for J = 0.93 (filled circles). For comparison 
rriT and 2m^ are also shown. 



23 



Gauss form, a slightly better discription is given by a peak with a finite width. We 
have therefore approximated the 5-function by a Breit-Wigner form 



- m L ) = - ■ {u _ mL)2 + {r/2)2 ■ (107) 
Is r default model function for MEM is then given by 

W=XiWw-m £ ). (108) 

For J = 0.92 and 0.93 we got the best results for a width T = 0.005, a maximal 
value uj m = 3, and a resolution Au = 0.001 of the frequency range. Since the Breit- 
Wigner form (107) is normalized to 1 in the range [— oo, oo] the cut-off at the limits 
of the u;-range [0, u m ] was corrected by a corresponding factor in A®. In Figs. 11 
and 12 we show the results for mi- The further approach of to the value of 2mr 
is manifest, in particular in Fig. 12. As expected from the result of Ref. [3], where 
on the critical line a ratio = 1.99(1) was obtained, we find at T c that 

coincides with Iwit- This is seen in Fig. 13, where we show the spectral functions 
Al(uj) (not Al{uj)) for each of our H- values. Again, we have used the default model 
of Eq. (108) with the width V = 0.005 and the resolution Au; = 0.001. The maximal 
value u m was in the range 2 — 4. 




Figure 13: The spectral function A l {oj) at T c and fixed values of H plotted versus 
H u °. For comparison the values of 2m-r (filled circles) are also shown. 
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One can use as well a larger width of for example 0.01 to obtain similarly good cor- 
relation functions. The peak heights shrink then correspondingly. It is remarkable 
that no further peak or continuum contributions appear in the spectral functions at 
the critical temperature. 

In the low temperature phase where we expect a continuum contribution in the 
longitudinal spectral function starting at the threshold value 2m-r we have tried 
several forms for the default model. First we have used the continuum form of Eq. 
(83). Here, the default spectral function is strictly zero below to = 2rriT and the 
MEM result remains always zero there. We have therefore moderated the model 
by combining the left half of a sharp Breit-Wigner form at the threshold with the 
continuum form above. That allows for non-zero spectral function values also below 
the threshold. Yet, this default model is predicting a too low threshold peak and 
too large contributions at high cu-values. It turned out that indeed a Breit-Wigner 
default model centered at the threshold nevertheless allows for continuum contri- 
butions well above its peak position and delivers the best results. In contrast, a 
corresponding Gaussian distribution leads to poor results. Obviously, the tail of the 
Gaussian distribution is dying out too fast, whereas the tail of the Breit-Wigner 
form falls off with u~ 2 just as the continuum default model. In Fig. 14 we show the 
MEM results for the spectral functions Al{uj) at J = 0.94, the temperature value 




Figure 14: The spectral function A L {u) at J = 0.94 and fixed values of H plotted 
versus if 1 / 2 . For comparison the values of 2m<r (filled circles) are also shown. 
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in the cold phase closest to T c . As in Fig. 6 we plot the data versus H l l 2 . The 
parameters we used for MEM and the default model were V = 0.01, Auj = 0.001 
and uj m was again in the range 2 — 4. We observe that apart from the main peak 
at the threshold value uj = 2m^ additional peaks appear at higher upvalues, in 
particular for small external fields H. In other respects the picture at J = 0.94 
resembles still the one at T c , Fig. 13: the peaks at the threshold are the major 
contributions to the spectral functions. 

The behaviour of the spectral function changes gradually when the temperature is 
shifted deeper into the cold phase. Apart from the smallest if -values the continuum 
contributions become more and more important. Correspondingly, the threshold 
peak decreases faster and/or merges with the next peak at higher u such that some- 
times no peak remains at the threshold. Also, the u-range where relevant continuum 
contributions appear increases with decreasing temperature. All this is demon- 
strated in Figs. 15, 16 and 17, where we show our results for J = 0.95, 0.97, 1.0 
and 1.2 with different cu-ranges from to 1-2. In the plots for J = 1.0 and 1.2 we 
have limited the peaks at H = 0.0001 to enlarge the other contributions for bet- 
ter visibility. The parameters for MEM and the default model were T = 0.01 and 
Alu = 0.001. Since uj m sometimes extended up to 6 we have also tried Au; = 0.0015 
or 0.002 and found no difference in the results. 



A L (u) 




Figure 15: The spectral function Al{uj) at J = 0.95 and fixed values of H plotted 
versus if 1//2 . For comparison the values of 2my (filled circles) are also shown. 
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Figure 17: The spectral function Al{uj) at J = 1.2 and fixed values of H plotted 
versus H 1 ' 2 . For comparison the values of 2mr (filled circles) are also shown. 



6 Additional results 
6.1 Test of the PP-relation 

In Section 5.1 we have shown that the transverse correlators have the Gaussian 
form and we have determined the respective transverse masses my. This enables us 
now to test the relation of Patashinskii and Pokrovskii which we have discussed in 
Section 2.4 . In Fig. 18 we plot the ratio of the measured longitudinal susceptibility 
Xl to the result Xl P °f Eq. (67) as a function of H for all our couplings J > J c 
in the low temperature phase. As expected, the ratio is close to 1 for very small 
if-values. For all couplings the ratio increases with increasing external field and 
reaches essentially a constant value. With decreasing temperature (increasing J) 
the ratio approaches 1 from above and at our lowest temperature T = 1/1.2 we 
are for all if- values already below 1.1 . The test of Xl P ls °f course a global one. 
More details are obtained from a direct comparison of the longitudinal correlator 
Dl(t) to its counterpart from Eq. (70). In Fig. 19 and 20 we show all our measured 
correlation functions for H = 0.0001 in the low temperature phase together with 
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Figure 18: The ratio Xl/Xl P °f the measured susceptibility xl to the result from 
the PP-relation, Eq. (67), for J = 0.94 (circles), 0.95 (diamonds), 0.97 (squares), 
1.0 (triangles) and 1.2 (davidstars) versus H. 
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Figure 19: The correlator Dl(t) for J = 0.94 (circles) and 0.95 (diamonds) at 
H = 0.0001. The lines show the respective results from the PP-relation (70). 
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Figure 20: The correlator D L {r) for J = 0.97 (squares), 1.0 (triangles) and 1.2 
(davidstars) at H = 0.0001. The lines show the respective results from the PP- 
relation (70). 

the respective results from (70). In Fig. 19 we plot the correlators for J = 0.94 
and 0.95, the two couplings closest to the critical point. Whereas for J = 0.94 the 
prediction from Eq. (70) is still definitely different from the data, the line and the 
data are close to each other already for J = 0.95, apart from the small r-region. In 
Fig. 20 we observe with decreasing temperature a further approach, at J = 1.2 the 
prediction coincides fully with the data. In the high temperature region the relation 
(50) obviously cannot hold for H — > : there is no spontaneous magnetization and 
the longitudinal and transverse correlation functions must agree with each other. 
Accordingly, xl an d Xl P can differ by a factor of 10 2 for small external fields. 

6.2 The stiffness 

In Sections 2.4 and 3 we have discussed the behaviour of the stiffness c s . From 
our results for the transverse mass ttit and the transverse susceptibility \t we can 
calculate c s using Eq. (51) . In Fig. 21 we show the stiffness as a function of H 1 / 2 (see 
Eq. (101) ) for fixed temperatures below T c . As expected, the stiffness increases with 
decreasing temperature. With decreasing external field the stiffness grows slightly, 
however for H —>■ a finite value will be reached. This is different on the critical 
line, where 

c s = H [gT(0)} 2 h~^ = 0.9740 H'^" , (109) 
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Figure 21: The stiffness c s versus if 1//2 in the low temperature phase for J = 0.94 
(circles), 0.95 (diamonds), 0.97 (squares), 1.0 (triangles) and 1.2 (davidstars). The 
lines are intended to guide the eye. 



and c s is weakly diverging for H — > 0. We have calculated the prefactor of H~ riVc 
from the critical parameters given in Ref. [3]. The prediction from the last equation 
is confirmed by our data as can be seen in Fig. 22. In the high temperature region 
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Figure 22: The stiffness c s on the critical line versus H nVc and the prediction (line) 
from Eq. (109). 
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Figure 23: The stiffness c s as a function of H in the high temperature phase for 
J = 0.93 (circles), 0.92 (diamonds), 0.91 (squares) and 0.90 (triangles). The dotted 
lines are drawn to guide the eye. 



c s is close to one and only slowly varying with the external field. This is evident 
from Fig. 23, where we have plotted function of H and not H 2 to avoid a 

too dense plot at small H. 



6.3 Scaling of the longitudinal correlation length 

In Ref. [3] the scaling function of the longitudinal correlation length was calculated 
from the directly measured values of For T >T C the longitudinal mass is the 
inverse of the longitudinal correlation length. Such an identification is not possible 
in the low temperature range where we have a whole continuum of contributing 
or ^-values. In the high temperature region however we can use the relation 

9i (z) = ^ , (110) 

to calculate the scaling function from our m^-data. The result is shown in Fig. 
24 together with the asymptotic form for large z- values from Eq. (96). Our data 
are scaling surprisingly well, better than the former results in Fig. 9 of Ref. [3], 
in particular in the small ^-region. One reason for the improvement may be the 
increase in the lattice size used, the other one is certainly the safer determination of 
rriL as compared to that of £5. As in Ref. [3] we note the similarity in the form of 
this scaling function and that of the longitudinal susceptibility. Consequently, both 
functions peak at about the same value of z, z v = 1.335 . 
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Figure 24: The normalized scaling function g^{z) in the high temperature phase for 
J = 0.93 (circles), 0.92 (diamonds), 0.91 (squares) and 0.90 (triangles). The dashed 
line indicates the peak position of xl, the solid line is the asymptotic form (96). 

6.4 Magnetization and susceptibilities at T c 

The behaviour of the magnetization M and the susceptibilities xl and xt on the 
critical line T = T c is controlled by the critical exponent 5 or rather its inverse as 
is evident from Eqs. (87) and (88). In Fig. 25 we show our data for M at T c as a 
function of if 1//<5 , where S = 4.824, the value of [3] was used. We compare the data 
with the result from Eq. (87) with B c = 0.721, again from [3]. Except for the point 
at H = 0.0001 which perhaps shows some finite size effect, all data agree very well 
with the line. A straight line fit to the data for In M as a function of InH without 
the first point leads to a somewhat more accurate result 

B c = 0.7198 ± 0.0005 , and 5 = 4.831 ± 0.003 . (Ill) 

The ratio Xt/Xl offers a straightforward way to determine 5: on the critical line 
it should coincide with the exponent. The drawback is however that the ratio has 
larger errors and that it shows definite finite size effects at small if-values. The 
latter are due to the longitudinal susceptibility. In Fig. 26 we show our data for the 
ratio. The values for H < 0.003 are still affected by the volume dependence. In 
order to determine S we have therefore used only the points from H > 0.0005. We 
find 

5 = 4.834 ±0.007 . (112) 
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Figure 25: The magnetization M (diamonds) as a function of if 1//<5 on the critical 
line. The dashed line is given by Eq. (87) with the parameters from Ref. [3]. 
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Figure 26: The ratio Xt/xl (diamonds) versus H on the critical line. The dashed 
line is 5 = 4.824 from Ref. [3], the solid line the fit result 5 = 4.834. 



7 Summary and conclusions 



In our paper we have investigated the transverse and longitudinal 2-plane spin cor- 
relation functions of the three-dimensional 0(4)-invariant non-linear a-model. The 
data were obtained from simulations on a cubic lattice with linear extension L = 120 
for small external fields H and temperatures below, above and at the critical tem- 
perature T c . The main objective was the determination of the spectral functions of 
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the correlators. In particular we were interested in the identification of the lowest 
states of the spectra tot,l, which relate to the correlation lengthes by = l/ m r,L, 
the influence of possible higher states on the correlation functions, and the interplay 
of the transverse and longitudinal channels. In order to calculate the spectral func- 
tions we used Bryan's algorithm for the maximum entropy method with a modified 
kernel which improves numerical stability at small u. 

As a result we find in the transverse case a single sharp peak in the spectral func- 
tion for each of our T and i/-values. The correlator is that of a free particle with 
mass tot, that is it has the known Gaussian form. We confirm therefore also corre- 
sponding assumptions of spin-wave theory about the dominance of long-wavelegth 
transverse fluctuations (i.e. small tot) of Gaussian type for small fields below T c . In 
the very high temperature region we observe a single sharp peak as well in the longi- 
tudinal spectrum. On approaching the critical point from above the peak broadens 
somewhat and at T c its position uil coincides with 2to^ for all our H— values. In 
the low temperature region T < T c we still find a significant peak at u = 2rriT- 
At higher u— values a continuum of states with several peaks of decreasing heights 
appears. This is expected from a relation of Patashinskii and Pokrovskii between 
the longitudinal and the transverse correlation functions. However, a comparison of 
the relation with the correlation function data reveals that there must be additional 
contributions to the longitudinal spectra at higher external fields. With decreasing 
temperatures these additive contributions gradually disappear. 

The influence of the longitudinal spin component on the transverse ones is charac- 
terized by the stiffness c s . As expected we find an increase of c s with decreasing 
temperature below T c . Above the critical temperature the stiffness is close to one 
and only weakly dependent on the external field. At T c we verify the predicted crit- 
ical behaviour c s ~ H~ VUc . With our results for to^ for the high temperature region 
we have calculated the scaling function of the longitudinal correlation length. The 
new data scale surprisingly well, better than those of Ref. [3], which were derived 
from direct measurements of the correlation length. As a last check we have com- 
pared the data for the magnetization and the susceptibilities to the predicted critical 
behaviour at T c . We find again agreement with the former values B c = 0.721(2) and 
5 = 4.824(9) from [3]. A new fit to our data for the magnetization leads however to 
the more accurate numbers B c = 0.7198(5) and 5 = 4.831(3). 
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